rm(list=ls())
setwd("/Users/yuehou/Dropbox/KingReplication/Ongoing Work")
 

#############################
#############################
#### Loop with Beck-Katz ####
#############################
#############################

# NOTE: You must input the models (below) manually 

# This is a hack of zelig functionality, that allows for correct
# standard errors and coefficients, but also allows the use of
# non-zelig models (plm), etc
# It also allows the user to modify the imputed dataset by adding lags, etc


load("imputations23.Rdata")
m <- 8

# Loads initial imputations into dataframe
for (j in 1:m){
	assign(paste("data_",j,sep=""), amelia.out9$imputation[[j]])
}


for (j in 1:m){
	data_temp <- get((paste("data_",j,sep="")))
	

	# Reframe as a panel dataset
	library(plm)
	p_newset <- pdata.frame(data_temp)
	
	p_newset2<-subset(data_temp,dhasminwage==1)  #just for minwage
	p_newset2<-pdata.frame(p_newset2)

	p_newset$log_oecd_gdp <- log(p_newset$oecd_gdp)

# Create lag terms here
	p_newset$lag_ia5010 <- lag(p_newset$final_ia5010)
	p_newset$lag_govem_yh2<- lag(p_newset$govem_yh2)

	p_newset$lag_leftc<- lag(p_newset$leftc)
	p_newset$lag_leftc_2<- lag(p_newset$leftc,2)
	p_newset$lag_leftc_3<- lag(p_newset$leftc,3)
	p_newset$lag_leftc_4<- lag(p_newset$leftc,4)
	p_newset$lag_leftc_5<- lag(p_newset$leftc,5)

	p_newset$lag_leftgs   <- lag(p_newset$leftgs)
	p_newset$lag_leftgs_2 <- lag(p_newset$leftgs,2)
	p_newset$lag_leftgs_3 <- lag(p_newset$leftgs,3)
	p_newset$lag_leftgs_4 <- lag(p_newset$leftgs,4)
	p_newset$lag_leftgs_5 <- lag(p_newset$leftgs,5)

	p_newset$lag_arm_generos = lag(p_newset$arm_generos)
	p_newset$lag_minwag_jf = lag(p_newset$minwag_jf)
	
	

# Transform generosity measure

	temp1 <- p_newset$soc_exp1 <- (p_newset$arm_generos - p_newset$lag_arm_generos)
	temp2 <- p_newset$lag_arm_generos
	temp3 <- p_newset$soc_exp <- 100 * (temp1 / temp2)
	
	p_newset$lag_soc_exp<- lag(p_newset$soc_exp)

	#assign(paste("as_",j,sep=""), p_newset$arm_generos)
	#assign(paste("as2_",j,sep=""), p_newset$lag_arm_generos)
	#assign(paste("as3_",j,sep=""), p_newset$soc_exp )


	# Oil covariate
	p_newset$oilshock <- 0
	p_newset$oilshock <- (factor(year)==1973 | factor(year) == 1972 | factor(year) == 1971 | 	factor(year) == 1970)



	
 ###############################################
   

    # import model
    
   #results<- plm(model="within",data=p_newset,final_ia5010 ~ minwag_jf + hkcorp_1:minwag_jf + dhasminwage + lag_ia5010 + govem_yh2 + govem_yh2:hkcorp_1 + soc_exp  + hkcorp_1:soc_exp + hkcorp_1 + arm_unemp + ldc_imports_scaled + oecd_gdpgr + dGermany + oilshock)
   
   results<-plm(model="within",data=p_newset,final_ia5010 ~ minwag_jf + vk_corp:minwag_jf + dhasminwage + lag_ia5010 + govem_yh2 + govem_yh2:vk_corp + soc_exp  + vk_corp:soc_exp + vk_corp + arm_unemp + ldc_imports_scaled + oecd_gdpgr + dGermany + oilshock)

   results2 <- coeftest(results, vcov=function(x) vcovBK(x, cluster=c("group","time")))
 
   # marginal effects
   hkvalue<-seq(1,5,1)
   
   beta_interaction<- results$coeff["govem_yh2"]+hkvalue*results$coeff["vk_corp:govem_yh2"]
   covariance1 <- vcovBK(results,cluster=c("group","time"))["govem_yh2","vk_corp:govem_yh2"]
   
      se_interaction <- sqrt(results2["govem_yh2",2]^2+(hkvalue^2)*results2["vk_corp:govem_yh2",2]^2+2*hkvalue*covariance1)
   
   
     
   results <- coeftest(results, vcov=function(x) vcovBK(x, cluster=c("group","time")))

	# Output vectors
	assign(paste("resultsc_",j,sep=""),results[,1])
	assign(paste("resultss_",j,sep=""),results[,2])
	assign(paste("resultsinteraction_",j,sep=""),beta_interaction)
	assign(paste("resultsinterSE_",j,sep=""),se_interaction)
}




# Vectors of Results 	
coef <- 0
coef_interaction<-0
se <- 0
se2 <-0
se_interaction<-0
se_interaction_2<-0

for (j in 1:m){
	coef <- coef + get(paste("resultsc_",j,sep=""))
	coef_interaction<- coef_interaction + get(paste("resultsinteraction_",j,sep=""))
	se <- se + get(paste("resultss_",j,sep="")) 
	se2 <- se2 + (get(paste("resultss_",j,sep="")))^2
	se_interaction<- se_interaction + get(paste("resultsinterSE_",j,sep="")) 
	
	se_interaction_2<- se_interaction_2 + (get(paste("resultsinterSE_",j,sep="")))^2
}


coef <- coef/m
coef_interaction<-coef_interaction/m
se <- se/m
se2 <- se2/m
se_interaction_2 <-se_interaction_2/m
se3 <- 0
se4 <- 0


for (j in 1:m){
	se3 <- se3 + (get(paste("resultsc_",j,sep="")) - coef)^2
	se4 <- se4 + (get(paste("resultsinteraction_",j,sep=""))-coef_interaction)^2
}


se3 <- se3/(m-1)
se4 <- se4/(m-1)
final_se <- sqrt(se2 + se3*(1+(1/m)))
t <- coef/final_se

final_se_interaction<-sqrt(se_interaction_2 + se4*(1+(1/m)))

t_interaction<-coef_interaction/final_se_interaction

# Output Results
cbind(coef,final_se,t)


cbind(coef_interaction,final_se_interaction,t_interaction)


####### PLOT
vkrange2 = seq(1,5,.01)
vkrange=c(1,2,3,4,5)
par(mfrow=c(1,3))
## iv = govem

plot(vkrange, coef_interaction,"l",lwd=2,ylim=c(-.1,.1),xlab="Corporatism",ylab="Marg. Effect of Public Sector",main="Inequality")
lines(lowess(vkrange,coef_interaction + 1.96*final_se_interaction),col="red",lwd=1,lty=3)
lines(lowess(vkrange,coef_interaction - 1.96*final_se_interaction),col="red",lwd=1,lty=3)
abline(h=0,col="blue",lwd=.5)





